A trivial SOS decomposition example

Contributed by: votroto Adapted from: Examples 3.25 of [BPT12]

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R. Semidefinite Optimization and Convex Algebraic Geometry. Society for Industrial and Applied Mathematics, 2012.

using DynamicPolynomials
using SumOfSquares
import CSDP

The polynomial p = x^2 - x*y^2 + y^4 + 1 is SOS. We can, for example, decompose it as p = 3/4*(x - y^2)^2 + 1/4*(x + y)^2 + 1, which clearly proves that p is SOS, and there are infinitely many other ways to decompose p into sums of squares.

We can use SumOfSquares.jl to find such decompositions.

First, setup the polynomial of interest.

@polyvar x y
p = x^2 - x*y^2 + y^4 + 1

\[ y^{4} - xy^{2} + x^{2} + 1 \]

Secondly, constrain the polynomial to be nonnegative. SumOfSquares.jl transparently reinterprets polyonmial nonnegativity as the appropriate SOS certificate for polynomials nonnegative on semialgebraic sets.

model = SOSModel(CSDP.Optimizer)
@constraint(model, cref, p >= 0)

cref : $ (1)y^{4} + (-1)xy^{2} + (1)x^{2} + (1) \text{ is SOS} $

Thirdly, optimize the feasibility problem!

optimize!(model)

Success: SDP solved
Primal objective value: 0.0000000e+00
Dual objective value: 0.0000000e+00
Relative primal infeasibility: 4.15e-17
Relative dual infeasibility: 5.00e-11
Real Relative Gap: 0.00e+00
XZ Relative Gap: 3.63e-10
DIMACS error measures: 6.34e-17 0.00e+00 5.00e-11 0.00e+00 0.00e+00 3.63e-10
CSDP 6.2.0
This is a pure primal feasibility problem.
Iter:  0 Ap: 0.00e+00 Pobj:  0.0000000e+00 Ad: 0.00e+00 Dobj:  0.0000000e+00
Iter:  1 Ap: 9.00e-01 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  2.7659722e+01
Iter:  2 Ap: 1.00e+00 Pobj:  0.0000000e+00 Ad: 1.00e+00 Dobj:  1.2562719e+01

Lastly, recover a SOS decomposition. In general, SOS decompositions are not unique!

sos_dec = sos_decomposition(cref, 1e-4)

(-0.9492757372229986*y^2 + 0.5916342145954683*x + 0.7423566403018638)^2 + (-1.1201589623631873*y)^2 + (5.434579573573638e-16*y^2 + 0.7820242435285943*x - 0.6232480104529254)^2 + (-0.31444486753600015*y^2 - 0.19597713808894568*x - 0.24590350966628965)^2

Converting, rounding, and simplifying - Huzza, Back where we began!

polynomial(sos_dec, Float32)

\[ y^{4} - xy^{2} + x^{2} - 4.440892e-16y^{2} + 1.9428903e-16x + 1.0 \]

A deeper explanation and the unexplained 1e-4 parameter

p = x^2 - x*y^2 + y^4 + 1 can be represented in terms of its Gram matrix as

gram = gram_matrix(cref)

GramMatrix{Float64, MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}, Float64, SymMatrix{Float64}}([1.0 -0.49999999999999994 0.0 -0.627378050481286; -0.49999999999999994 1.0 0.0 2.331468351712827e-17; 0.0 0.0 1.2547561009625725 0.0; -0.627378050481286 2.331468351712827e-17 0.0 1.0], MonomialBasis{DynamicPolynomials.Monomial{true}, DynamicPolynomials.MonomialVector{true}}(DynamicPolynomials.Monomial{true}[y², x, y, 1]))

gram.basis.monomials' * gram.Q * gram.basis.monomials

\[ y^{4} - 0.9999999999999999xy^{2} + x^{2} + 4.440892098500626e-16y^{2} + 4.662936703425654e-17x + 1.0 \]

where the matrix gram.Q is positive semidefinite, because p is SOS. If we could only get the decomposition gram.Q = V' * V, the SOS decomposition would simply be ||V * monomials||^2.

Unfortunately, we can not use Cholesky decomposition, since gram Q is only semidefinite, not definite. Hence, SumOfSquares.jl uses SVD decomposition instead and discards small singular values (in our case 1e-4).

This page was generated using Literate.jl.